* This Stata dofile is one of two that is written to accompany the paper:
* J Gans & A Leigh, 'Born on the First of July: An (Un)natural Experiment in Birth Timing' (2009) Journal of Public Economics, 93(1-2): 246-263 
* Feel free to use or adapt it, but please cite that paper.
* Questions to andrew_leigh@ksg02.harvard.edu

* Note that we are not permitted to share the data from this paper. Other users should contact the 
* Australian Bureau of Statistics and the Australian Institute of Health and Welfare for the respective datasets.
* Note that there are vintage effects in the ABS data, and our ABS births data are births registered by 31 Dec 2006.

version 10

*******************************
* Baby Bonus as a percent of AWE
* Using 2004 HILDA data.
*******************************
use dhifdip dhifdin dhhwte using "C:\Users\Andrew\Datasets\HILDA\combined_d40c.dta", clear
gen hhdispinc=dhifdip-dhifdin
sum hhdispinc [aw=dhhwte],d
di (3000/r(p50))
di 3000/55511
di (3000/55511)*52

use ehifdip ehifdin ehhwte using "C:\Users\Andrew\Datasets\HILDA\Combined_e51c.dta", clear
gen hhdispinc=(ehifdip-ehifdin)*1.04
sum hhdispinc [aw=ehhwte],d
di (834/r(p50))
di (834/r(p50))*52
di (834/r(p25))*52

*******************************
* Setting up AIHW data
*******************************
cd "C:\Users\Andrew\My publications\Aust - baby bonus bunching\"
use births_aihw.dta, clear
for any "Vaginal, not induced" "Vaginal,induced" "c-section, induced" "c-section, not induced" \ any vni vi csi csni: replace birth_type="Y" if birth_type=="X"
reshape wide seps, i(date) j(birth_type) str
renpfix seps
sort date
save temp_births_aihw, replace

*******************************
* ABS regressions (incl 2006) 
*******************************
clear
set more off
set mem 50m
set matsize 800
cd "C:\Users\Andrew\My publications\Aust - baby bonus bunching\"
use births_by_day_2006.dta,clear
renpfix _ births
reshape long births,i(month day) j(year)
drop if day==0
ren month monthname
gen month=.
for num 1/12 \ any January February March April May June July August September October November December: replace month=X if monthname=="Y"
gen date=mdy(month,day,year)
format date %d
gen dow=dow(date)
gen halfyear=halfyear(date)
* The day of year shouldn't move in leap years
gen doy=doy(date)
replace doy=doy+1 if year/4~=int(year/4) & month>2
sum doy if month==7 & day==1
gen babybonus=1 if halfyear==2 & year==2004
recode babybonus .=0
gen babybonus2=1 if halfyear==2 & year==2006
recode babybonus2 .=0
gen logbirths=log(births)
gen qb=1 if dow==1 & month==6 & day>=8 & day<=14
recode qb .=0

* Dropping data after Aug 2006 (since reports start to decline)
drop if month>8 & year==2006

* Merging on the AIHW births data
sort date
merge date using temp_births_aihw
tab _merge 
list date if _merge==2
drop _merge
gen births_aihw=vni+vi+csi+csni
gen births_nonnatural=(vi+csi+csni)/(vni+vi+csi+csni)
list date births births_aihw if births~=. & births_aihw~=.

* Where are 29-30 June and 1-2 July in the data?
drop if date==.
xtile births_p=births, nq(100)
xtile births_p_nw=births if dow>=1 & dow<=5, nq(100)
xtile births_p_w=births if dow==0 | dow==6, nq(100)
sort doy
list day month doy births births_p births_p_nw births_aihw if doy>=180 & doy<=185 & year==2004
list day month doy births births_p births_p_w births_p_nw if doy>=180 & doy<=185 & year==2006
sum births if dow>=1 & dow<=5
sum births if dow==0 | dow==6

gen holiday=0
replace holiday=1 if (day==1 & month==1) | (day==25 & month==12) | (day==26 & month==12) | (day==26 & month==1) | (day==25 & month==4)
sort date
merge date using easter_dates.dta, nokeep
drop _merge
tsset date
replace holiday=1 if f2.easter==1 | l.easter==1
replace holiday=1 if qb==1

gen births_p_d=.
for num 0/6: xtile temp=births if year==2006 & dow==X, nq(100) \ replace births_p_d=temp if dow==X \ drop temp
list day month doy births births_p_d if doy>=175 & doy<=190 & year==2006

* Ranking births
bysort dow: egen birthsrank=rank(births), field
list day month doy births birthsrank if doy>=175 & doy<=190 & year==2006

* Comparing weekly results
for num 2004/2006: egen birthsweek1=sum(births) if month==6 & day>=24 & year==X \ egen birthsweek2=sum(births) if month==7 & day<=7 & year==X \ sum birthsweek* \ drop birthsweek*

gen window7=1 if (month==6 & day>=24) | (month==7 & day<=7)
gen window14=1 if (month==6 & day>=17) | (month==7 & day<=14)
gen window21=1 if (month==6 & day>=10) | (month==7 & day<=21)
gen window28=1 if (month==6 & day>=3) | (month==7 & day<=28)
gen window56=1 if doy>=182.5-56 & doy<=182.5+56

* Main results
for num 7: xi: areg births babybonus2 qb i.year*i.dow if windowX==1 & year~=2004, a(doy) \ outreg using results_babybonus2006.doc, coefastr nocons bracket 3aster replace bdec(3) se ct("Births - X day window") addstat("Moved 2006",_b[babybonus2]*X/2) adec(0) 
for num 14 21 28 : xi: areg births babybonus2 qb i.year*i.dow if windowX==1 & year~=2004, a(doy) \ outreg using results_babybonus2006.doc, coefastr nocons bracket 3aster append bdec(3) se ct("Births - X day window")  addstat("Moved 2006",_b[babybonus2]*X/2) adec(0)
for num 7 14 21 28 : xi: areg logbirths babybonus2 qb i.year*i.dow if windowX==1 & year~=2004, a(doy) \ outreg using results_babybonus2006.doc, coefastr nocons bracket 3aster append bdec(3) se ct("Ln(Births) - X day window") addstat("Moved 2006",100*(exp(_b[babybonus2]/2)-1)) adec(0)

* Modelling all the small changes too
gen bb_amount=0 
for any 1jul2004 20sep2004 20mar2005 20sep2005 20mar2006 01jul2006 20sep2006 20mar2007 \ num 3000 3042 3079 3119 3166 4000 4100 4133: replace bb_amount=Y if date>=d(X) & date~=.
for any 20sep2004 20mar2005 20sep2005 20mar2006 20sep2006 20mar2007 \ any sep04 mar05 sep05 mar06 sep06 mar07: gen bb_Y=0 \ replace bb_Y=1 if date>=d(X) & date<=d(X)+28 
replace window7=1 if (month==3 & day>=13 & day<=26) | (month==9 & day>=13 & day<=26)
xi: areg births babybonus babybonus2 bb_sep04 bb_mar05 bb_sep05 bb_mar06 qb i.year*i.dow if window7==1, a(doy)
xi: areg births bb_sep04 bb_mar05 bb_sep05 bb_mar06 i.year*i.dow if window7==1 & (month==3 | month==9), a(doy)
xi: areg logbirths babybonus babybonus2 bb_sep04 bb_mar05 bb_sep05 bb_mar06 qb i.year*i.dow if window7==1, a(doy)
xi: areg logbirths bb_sep04 bb_mar05 bb_sep05 bb_mar06 i.year*i.dow if window7==1 & (month==3 | month==9), a(doy)
xi: areg logbirths bb_amount i.year*i.dow if window7==1 & (month==3 | month==9), a(doy)

* Graphing the ABS series
bysort doy: egen births_mean=mean(births) if year<=2003
set scheme s2mono
gen bb_period=0
replace bb_period=1 if (month==6 & year==2004) | (month==7 & year==2004) | (month==6 & year==2006) | (month==7 & year==2006)
xi: reg births i.dow*i.year holiday if bb_period==0
predict births_dow
gen births_sa=births-births_dow
gen births_sa2=births_sa if doy<=182 & births_sa<0
replace births_sa2=births_sa if doy>=183 & births_sa>0

tabstat births_dow,by(dow)
twoway line births doy if doy>=155 & doy<=210 & year==2006,yline(0) legend(off) xline(182.5) xtitle("") yti("Births per day (unadjusted)") xlabel(155 "June 3" 162 "June 10" 169 "June 17" 176 "June 24" 183 "July 1" 189 "July 7" 196 "July 14" 203 "July 21" 210 "July 28",angle(vertical)) ti("Panel A: Raw Data") name(raw,replace) nodraw
twoway area births_sa2 doy if doy>=155 & doy<=210 & year==2006, color(gray) base(0) cmiss(n) || line births_sa doy if doy>=155 & doy<=210 & year==2006,yline(0) legend(off) xline(182.5) xtitle("") yti("Births per day (relative to expected)") xlabel(155 "June 3" 162 "June 10" 169 "June 17" 176 "June 24" 183 "July 1" 189 "July 7" 196 "July 14" 203 "July 21" 210 "July 28",angle(vertical)) ti("Panel B: Regression-Adjusted") name(adj,replace) nodraw
gr combine raw adj, ti("Figure 8: The Introduction Effect in 2006") cols(1) altshrink copies note("Panel A shows daily birth counts. Panel B shows births relative to expected, accounting for" "year, day of week, day of year, and public holidays." "In Panel B, shading shows days of unusually low births in June and unusually high births in July.")

twoway line births doy if doy>=155 & doy<=210 & year==2005,yline(0) legend(off) xline(182.5) xtitle("") yti("Births per day (unadjusted)") xlabel(155 "June 3" 162 "June 10" 169 "June 17" 176 "June 24" 183 "July 1" 189 "July 7" 196 "July 14" 203 "July 21" 210 "July 28",angle(vertical)) ti("Panel A: Raw Data") name(raw,replace) nodraw
twoway area births_sa2 doy if doy>=155 & doy<=210 & year==2005, color(gray) base(0) cmiss(n) || line births_sa doy if doy>=155 & doy<=210 & year==2005,yline(0) legend(off) xline(182.5) xtitle("") yti("Births per day (relative to expected)") xlabel(155 "June 3" 162 "June 10" 169 "June 17" 176 "June 24" 183 "July 1" 189 "July 7" 196 "July 14" 203 "July 21" 210 "July 28",angle(vertical)) ti("Panel B: Regression-Adjusted") name(adj,replace) nodraw
gr combine raw adj, ti("Appendix Figure 1: The Introduction Effect in 2005") subti("(Falsification Exercise)") cols(1) altshrink copies note("Panel A shows daily birth counts. Panel B shows births relative to expected, accounting for" "year, day of week, day of year, and public holidays." "In Panel B, shading shows days of unusually low births in June and unusually high births in July.")

twoway area births_sa2 doy if doy>=155 & doy<=210 & year==2004, color(gray) base(0) cmiss(n) || line births_sa doy if doy>=155 & doy<=210 & year==2004,yline(0) legend(off) xline(182.5) xtitle("") yti("Births per day (relative to expected)") xlabel(155 "June 3" 162 "June 10" 169 "June 17" 176 "June 24" 183 "July 1" 189 "July 7" 196 "July 14" 203 "July 21" 210 "July 28",angle(vertical)) ti("2004") name(adj1,replace) nodraw
twoway area births_sa2 doy if doy>=155 & doy<=210 & year==2005, color(gray) base(0) cmiss(n) || line births_sa doy if doy>=155 & doy<=210 & year==2005,yline(0) legend(off) xline(182.5) xtitle("") yti("Births per day (relative to expected)") xlabel(155 "June 3" 162 "June 10" 169 "June 17" 176 "June 24" 183 "July 1" 189 "July 7" 196 "July 14" 203 "July 21" 210 "July 28",angle(vertical)) ti("2005") name(adj2,replace) nodraw
twoway area births_sa2 doy if doy>=155 & doy<=210 & year==2006, color(gray) base(0) cmiss(n) || line births_sa doy if doy>=155 & doy<=210 & year==2006,yline(0) legend(off) xline(182.5) xtitle("") yti("Births per day (relative to expected)") xlabel(155 "June 3" 162 "June 10" 169 "June 17" 176 "June 24" 183 "July 1" 189 "July 7" 196 "July 14" 203 "July 21" 210 "July 28",angle(vertical)) ti("2006") name(adj3,replace) nodraw
gr combine adj1 adj2 adj3, ti("Appendix Figure 1: Including 2005 as a Placebo") altshrink copies note("Shading shows days of unusually low births in June and unusually high births in July")

*******************************
* ABS regressions (up to 2004)
*******************************
clear
set more off
set mem 150m
set matsize 800
cd "C:\Users\Andrew\My publications\Aust - baby bonus bunching\"
use births_by_day_2006.dta,clear
renpfix _ births
reshape long births,i(month day) j(year)
ren month monthname
gen month=.
for num 1/12 \ any January February March April May June July August September October November December: replace month=X if monthname=="Y"
gen date=mdy(month,day,year)
format date %d

* Drop 2005 & 2006
drop if year==2005 | year==2006
bysort year: egen annualbirths=sum(births)
tabstat annualbirths, by(year)
drop annualbirths

gen dow=dow(date)
gen halfyear=halfyear(date)
* The day of year shouldn't move in leap years
gen doy=doy(date)
replace doy=doy+1 if year/4~=int(year/4) & month>2
sum doy if month==7 & day==1
gen babybonus=1 if halfyear==2 & year==2004
recode babybonus .=0
gen logbirths=log(births)
gen qb=1 if dow==1 & month==6 & day>=8 & day<=14
recode qb .=0

* Merging on the AIHW births data
sort date
merge date using temp_births_aihw
tab _merge 
list date if _merge==2
drop _merge
gen births_aihw=vni+vi+csi+csni
gen births_nonnatural=(vi+csi+csni)/(vni+vi+csi+csni)
list date births births_aihw if births~=. & births_aihw~=.

* Where are 29-30 June and 1-2 July in the data?
drop if date==.
xtile births_p=births, nq(100)
xtile births_p_nw=births if dow>=1 & dow<=5, nq(100)
sort doy
egen births_rank=rank(births), field
list day month doy births births_p births_p_nw births_rank if doy>=180 & doy<=185 & year==2004
sum births if dow>=1 & dow<=5
sum births if dow==0 | dow==6

gen holiday=0
replace holiday=1 if (day==1 & month==1) | (day==25 & month==12) | (day==26 & month==12) | (day==26 & month==1) | (day==25 & month==4)
sort date
merge date using easter_dates.dta, nokeep
drop _merge
tsset date
replace holiday=1 if f2.easter==1 | l.easter==1
replace holiday=1 if qb==1

gen window7=1 if (month==6 & day>=24) | (month==7 & day<=7)
gen window14=1 if (month==6 & day>=17) | (month==7 & day<=14)
gen window21=1 if (month==6 & day>=10) | (month==7 & day<=21)
gen window28=1 if (month==6 & day>=3) | (month==7 & day<=28)
for any 35 56 112: gen windowX=1 if doy>=182.5-X & doy<=182.5+X

* Main results
for num 7: xi: areg births babybonus qb i.year*i.dow if windowX==1, a(doy) \ outreg using results_babybonus.doc, coefastr nocons bracket 3aster replace bdec(3) se ct("Births - X day window") addstat("Moved",_b[babybonus]*X/2) adec(0) 
for num 14 21 28 : xi: areg births babybonus qb i.year*i.dow if windowX==1, a(doy) \ outreg using results_babybonus.doc, coefastr nocons bracket 3aster append bdec(3) se ct("Births - X day window")  addstat("Moved",_b[babybonus]*X/2) adec(0)
for num 7 14 21 28 : xi: areg logbirths babybonus qb i.year*i.dow if windowX==1, a(doy) \ outreg using results_babybonus.doc, coefastr nocons bracket 3aster append bdec(3) se ct("Ln(Births) - X day window") addstat("Moved",100*(exp(_b[babybonus]/2)-1)) adec(0)

* Daily dummies
for num 7 14 21 28: gen beforeX=0 \ gen afterX=0
for num 155/182 \ num 28 28 28 28 28 28 28 21 21 21 21 21 21 21 14 14 14 14 14 14 14 7 7 7 7 7 7 7: replace beforeY=1 if year==2004 & doy==X
for num 183/210 \ num 7 7 7 7 7 7 7 14 14 14 14 14 14 14 21 21 21 21 21 21 21 28 28 28 28 28 28 28: replace afterY=1 if year==2004 & doy==X

for any births : xi: areg X before28 before21 before14 before7 after7 after14 after21 after28 holiday i.year*i.dow if window56==1,a(doy) \ outreg using results_babybonus.doc, coefastr nocons bracket 3aster append bdec(3) se ct("X - nonlinear")
for any before after: lincom 7*(X28+X21+X14+X7)
for any before after: lincom 7*(X21+X14+X7)
for any before after: lincom 7*(X14+X7)
for any before after: lincom 7*(X7)
for any logbirths: xi: areg X before28 before21 before14 before7 after7 after14 after21 after28 holiday i.year*i.dow if window56==1,a(doy) \ outreg using results_babybonus.doc, coefastr nocons bracket 3aster append bdec(3) se ct("X - nonlinear")

* Robustness check - 112 day window?
for any births : xi: areg X before28 before21 before14 before7 after7 after14 after21 after28 holiday i.year*i.dow if window112==1,a(doy) \ outreg using results_babybonus.doc, coefastr nocons bracket 3aster append bdec(3) se ct("X - nonlinear")
for any before after: lincom 7*(X28+X21+X14+X7)
for any before after: lincom 7*(X21+X14+X7)
for any before after: lincom 7*(X14+X7)
for any before after: lincom 7*(X7)

* Graphing the ABS series
set scheme s2mono
gen bb_period=0
replace bb_period=1 if (month==6 & year==2004) | (month==7 & year==2004) 
xi: reg births i.dow*i.year i.doy holiday if bb_period==0
predict births_dow
gen births_sa=births-births_dow
gen births_sa2=births_sa if doy<=182 & births_sa<0
replace births_sa2=births_sa if doy>=183 & births_sa>0
twoway line births doy if doy>=155 & doy<=210 & year==2004,yline(0) legend(off) xline(182.5) xtitle("") yti("Births per day (unadjusted)") xlabel(155 "June 3" 162 "June 10" 169 "June 17" 176 "June 24" 183 "July 1" 189 "July 7" 196 "July 14" 203 "July 21" 210 "July 28",angle(vertical)) ti("Panel A: Raw Data") name(raw,replace) nodraw
twoway area births_sa2 doy if doy>=155 & doy<=210 & year==2004, color(gray) base(0) cmiss(n) || line births_sa doy if doy>=155 & doy<=210 & year==2004,yline(0) legend(off) xline(182.5) xtitle("") yti("Births per day (relative to expected)") xlabel(155 "June 3" 162 "June 10" 169 "June 17" 176 "June 24" 183 "July 1" 189 "July 7" 196 "July 14" 203 "July 21" 210 "July 28",angle(vertical)) ti("Panel B: Regression-Adjusted") name(adj,replace) nodraw
gr combine raw adj, ti("Figure 1: The Introduction Effect in 2004") cols(1) altshrink copies note("Panel A shows daily birth counts. Panel B shows births relative to expected, accounting for" "year, day of week, day of year, and public holidays." "In Panel B, shading shows days of unusually low births in June and unusually high births in July.")
drop births_dow births_sa*

*****************************
* AIHW data
*****************************
* How closely do the series track one another?
set scheme s2mono
gen n1=_n if _n<1100 & _n>350
gen n2=_n if _n<1100 & _n>350
gen ax1=950 if _n==1
gen ay1=1067 if _n==1
gen ax2=990 if _n==1
gen ay2=1067 if _n==1
tw scatter births_aihw births if year==2004,msize(small) || line n1 n2, lpattern(solid) || pcarrow ay1 ax1 ay2 ax2, xlab(400(200)1000) ylab(400(200)1000) xti("ABS series (births registries)") yti("AIHW series (hospital records)") ti("Figure 2: Comparing Two Sources of Births Data") ysc(r(350 1100)) xsc(r(350 1100)) legend(off) note("Note: Dots show daily birth counts for 2004. The solid line is at x=y") text(1067 940 "July 1, 2004", place(w))
drop n1 n2 ay* ax*
* Regressing the ABS series on the AIHW series
reg births births_aihw if year==2004

* Share of each type of birth
for any vni vi csi csni: egen tempX=sum(X) if year==2004 \ egen temp_tot=sum(births_aihw) if year==2004 \ gen s_X=tempX/temp_tot \ drop temp*
sum s_*

* Graphing AIHW data
#delimit ;
for X in any vni vi csni csi \ Y in any "Vaginal, Not Induced" "Vaginal, Induced" "Cesarean Section, Not Induced" "Cesarean Section, Induced" \ Z in any A B C D:
gen births_X=X \ 
xi: reg births_X i.dow holiday if bb_period==0 \ 
predict births_dow \
gen births_sa=births_X-births_dow \
gen births_sa2=births_sa if doy<=182 & births_sa<0 \
replace births_sa2=births_sa if doy>=183 & births_sa>0 \
twoway area births_sa2 doy if doy>=155 & doy<=210 & year==2004, color(gray) base(0) cmiss(n) || line births_sa doy if doy>=155 & doy<=210 & year==2004,yline(0) legend(off) xline(182.5) xtitle("") yti("Births per day (relative to expected)") xlabel(155 "June 3" 162 "June 10" 169 "June 17" 176 "June 24" 183 "July 1" 189 "July 7" 196 "July 14" 203 "July 21" 210 "July 28",angle(vertical)) ti("Panel Z: Y") name(X,replace) nodraw \
drop births_dow births_sa*;
#delimit cr
set scheme s2mono
* Fig 3
twoway line births_vni doy if doy>=155 & doy<=210 & year==2004,legend(lab(1 "Vaginal, Not Induced")) || line births_vi doy if doy>=155 & doy<=210 & year==2004,legend(lab(2 "Vaginal, Induced")) || line births_csni doy if doy>=155 & doy<=210 & year==2004, legend(lab(3 "Cesarean Section, Not Induced")) || line births_csi doy if doy>=155 & doy<=210 & year==2004, legend(lab(4 "Cesarean Section, Induced")) legend(col(1)) xline(182.5) xtitle("") yti("Births per day (unadjusted)") xlabel(155 "June 3" 162 "June 10" 169 "June 17" 176 "June 24" 183 "July 1" 189 "July 7" 196 "July 14" 203 "July 21" 210 "July 28",angle(vertical)) ti("Figure 3: Comparing Birth Procedures") subti("Raw Data") note("Data cover June-July 2004, and show daily birth counts.") 
* Fig 4
gr combine vni vi csni csi, ti("Figure 4: Comparing Birth Procedures") subti("Regression-adjusted") altshrink note("Data cover June-July 2004, and show births relative to expected, accounting for day of week and public holidays." "Shading shows days of unusually low births in June and unusually high births in July.")

* Checking effects by birth procedure (vni vi csi csni)
for any vni vi csi csni: replace X=. if year~=2004
gen b=vni 
for num 7: xi: reg b babybonus holiday i.dow if windowX==1,  \ outreg using results_babybonus_aihw.doc, coefastr nocons bracket 3aster replace bdec(3) se ct("vni - X day window") addstat("Moved",_b[babybonus]*X/2) adec(0)
for num 14 21 28: xi: reg b babybonus holiday i.dow if windowX==1, \ outreg using results_babybonus_aihw.doc, coefastr nocons bracket 3aster append bdec(3) se ct("vni - X day window") addstat("Moved",_b[babybonus]*X/2) adec(0)
replace b=vi
for num 7 14 21 28: xi: reg b babybonus holiday i.dow if windowX==1, \ outreg using results_babybonus_aihw.doc, coefastr nocons bracket 3aster append bdec(3) se ct("vi - X day window") addstat("Moved",_b[babybonus]*X/2) adec(0)
replace b=csni 
for num 7 14 21 28: xi: reg b babybonus holiday i.dow if windowX==1, \ outreg using results_babybonus_aihw.doc, coefastr nocons bracket 3aster append bdec(3) se ct("csni - X day window") addstat("Moved",_b[babybonus]*X/2) adec(0)
replace b=csi 
for num 7 14 21 28: xi: reg b babybonus holiday i.dow if windowX==1, \ outreg using results_babybonus_aihw.doc, coefastr nocons bracket 3aster append bdec(3) se ct("csi - X day window") addstat("Moved",_b[babybonus]*X/2) adec(0) 
replace b=ln(vni) 
for num 7 14 21 28: xi: reg b babybonus holiday i.dow if windowX==1, \ outreg using results_babybonus_aihw.doc, coefastr nocons bracket 3aster append bdec(3) se ct("Ln(vni) - X day window") addstat("Moved",100*(exp(_b[babybonus]/2)-1)) adec(0)
replace b=ln(vi)
for num 7 14 21 28: xi: reg b babybonus holiday i.dow if windowX==1, \ outreg using results_babybonus_aihw.doc, coefastr nocons bracket 3aster append bdec(3) se ct("Ln(vi) - X day window") addstat("Moved",100*(exp(_b[babybonus]/2)-1)) adec(0)
replace b=ln(csni)
for num 7 14 21 28: xi: reg b babybonus holiday i.dow if windowX==1, \ outreg using results_babybonus_aihw.doc, coefastr nocons bracket 3aster append bdec(3) se ct("Ln(csni) - X day window") addstat("Moved",100*(exp(_b[babybonus]/2)-1)) adec(0)
replace b=ln(csi)
for num 7 14 21 28: xi: reg b babybonus holiday i.dow if windowX==1, \ outreg using results_babybonus_aihw.doc, coefastr nocons bracket 3aster append bdec(3) se ct("Ln(csi) - X day window") addstat("Moved",100*(exp(_b[babybonus]/2)-1)) adec(0)

replace b=vni 
for num 56: qui xi: reg b before28 before21 before14 before7 after7 after14 after21 after28 holiday i.dow if windowX==1, \ outreg using results_babybonus_aihw.doc, coefastr nocons bracket 3aster append bdec(3) se ct("vni - X day window") 
for any before after: lincom 7*(X28+X21+X14+X7)
replace b=vi
for num 56: qui xi: reg b before28 before21 before14 before7 after7 after14 after21 after28 holiday i.dow if windowX==1, \ outreg using results_babybonus_aihw.doc, coefastr nocons bracket 3aster append bdec(3) se ct("vi - X day window") 
for any before after: lincom 7*(X28+X21+X14+X7)
replace b=csni 
for num 56: qui xi: reg b before28 before21 before14 before7 after7 after14 after21 after28 holiday i.dow if windowX==1, \ outreg using results_babybonus_aihw.doc, coefastr nocons bracket 3aster append bdec(3) se ct("csni - X day window") 
for any before after: lincom 7*(X28+X21+X14+X7)
replace b=csi 
for num 56: qui xi: reg b before28 before21 before14 before7 after7 after14 after21 after28 holiday i.dow if windowX==1, \ outreg using results_babybonus_aihw.doc, coefastr nocons bracket 3aster append bdec(3) se ct("csi - X day window") 
for any before after: lincom 7*(X28+X21+X14+X7)
replace b=ln(vni)
for num 56: qui xi: reg b before28 before21 before14 before7 after7 after14 after21 after28 holiday i.dow if windowX==1, \ outreg using results_babybonus_aihw.doc, coefastr nocons bracket 3aster append bdec(3) se ct("Ln vni - X day window") 
replace b=ln(vi)
for num 56: qui xi: reg b before28 before21 before14 before7 after7 after14 after21 after28 holiday i.dow if windowX==1, \ outreg using results_babybonus_aihw.doc, coefastr nocons bracket 3aster append bdec(3) se ct("Ln vi - X day window") 
replace b=ln(csni) 
for num 56: qui xi: reg b before28 before21 before14 before7 after7 after14 after21 after28 holiday i.dow if windowX==1, \ outreg using results_babybonus_aihw.doc, coefastr nocons bracket 3aster append bdec(3) se ct("Ln csni - X day window") 
replace b=ln(csi) 
for num 56: qui xi: reg b before28 before21 before14 before7 after7 after14 after21 after28 holiday i.dow if windowX==1, \ outreg using results_babybonus_aihw.doc, coefastr nocons bracket 3aster append bdec(3) se ct("Ln csi - X day window") 

* Graph - share of Births That Are Non-Natural
gen births_nn=(vi+csi+csni)/(vni+vi+csi+csni)
xi: reg births_nn i.dow holiday if bb_period==0
predict births_nn_dow
gen births_nn_sa=births_nn-births_nn_dow
twoway line births_nn doy if doy>=155 & doy<=210 & year==2004, lpattern(solid) xline(182.5) xtitle("") yti("Share (unadjusted)") xlabel(155 "June 3" 162 "June 10" 169 "June 17" 176 "June 24" 183 "July 1" 189 "July 7" 196 "July 14" 203 "July 21" 210 "July 28",angle(vertical)) ti("Panel A: Raw Data") legend(off) name(raw,replace) nodraw
twoway line births_nn_sa doy if doy>=155 & doy<=210 & year==2004, lpattern(solid) yline(0) xline(182.5) xtitle("") yti("Share (relative to expected)") xlabel(155 "June 3" 162 "June 10" 169 "June 17" 176 "June 24" 183 "July 1" 189 "July 7" 196 "July 14" 203 "July 21" 210 "July 28",angle(vertical)) ti("Panel B: Regression-adjusted") legend(off) name(adj,replace) nodraw
gr combine raw adj, ti("Figure 5: Share of Births Delivered" "by Inducement and/or Cesarean Section") cols(1) altshrink copies note("Data cover June-July 2004. Panel A shows the unadjusted share, while Panel B shows the share" "relative to expected, accounting for day of week and public holidays.") 
tabstat births_nn if window7==1,by(month)
tabstat births_nn if window28==1,by(month)

